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Quantum Monte Carlo impurity solver to study the real-time dynamics of a Hubbard model which 
is driven out of equilibrium by a sudden increase in the on-site repulsion U. We discuss the im- 
plementation of the self-consistency procedure and some important technical improvements of the 
QMC method. The exact numerical solution is compared to iterated perturbation theory, which is 
found to produce accurate results only for weak interaction or short times. Furthermore we calculate 
the spectral functions and the optical conductivity from a Fourier transform on the finite Keldysh 
contour, for which the numerically accessible timescales allow to resolve the formation of Hubbard 
bands and a gap in the strongly interacting regime. The spectral function, and all one-particle 
quantities that can be calculated from it, thermalize rapidly at the transition between qualitatively 
different weak- and strong-coupling relaxation regimes. 
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I. INTRODUCTION 

The recent realization of a Mott insulating state of 
repulsively interacting fermions in trapped ultracold 
atomai*^ opens the door to controlled studies of the non- 
equilibrium properties of fermionic lattice models. At 
the same time, the relaxation dynamics of strongly cor- 
related electron systems is starting to be explored exper- 
imentally through femtosecond spectroscopy.—'^'^ Dy- 
namical mean-field theory! (DMFT) is a promising tool 
to approach these challenging issues from the theoretical 
side. The DMFT formalism is based on the mapping of 
a lattice model to a quantum impurity model. This ap- 
proximation is based on a purely spatial argument which 
becomes exact in the limit of infinite dimensions^ On 
the one hand this fact makes DMFT a nonperturbative 
method which can capture, e.g., the local Mott physics of 
the Hubbard model. On the other hand, it implies that 
DMFT can be formulated equally well in imaginary and 
real time, and hence the method can be applied to both 
equilibrium and nonequilibrium situationSjS 

A number of authors have employed the non- 
equilibrium DMFT framework to study dynamical prop- 
erties of the Falicov-Kimball model, which is a variant 
of the Hubbard model in which only one spin species 
can hop between lattice sites. Despite this simplifi- 
cation, the Falicov-Kimball model exhibits a relatively 
rich phase diagram with metallic, Mott-insulating, and 
charge-ordered phasesJ^ Its most attractive feature in 
the present context is that the associated quantum im- 
purity model in DMFT can be solved analytically or nu- 
merically by simple means which provides reliable ac- 
cess to the long-time dynamics. Both the transient dy- 
namics after the sudden switching-on of a static electric 
field^'^'^ and, using a combined Floquet and DMFT 
formalism, the non-equilibrium steady state in the pres- 



ence of an alternating or constant fieldi^'i^'ii were calcu- 
lated. The evolution of the momentum distribution and 
double occupation after an interaction quench, i.e., a sud- 
den change in the interaction parameter, was studied in 
Rcf. |T^, where it was also shown that for the Falicov- 
Kimball model these quantities do not thermalize. This 
is a consequence of the immobility of one spin species and 
the resulting quadratic form of the Hamiltonian for the 
mobile other spin species. 

A more realistic model for the description of correlated 
electron systems and interacting fermions in optical lat- 
tices is the Hubbard model. 



+ U{t)J2{nrT-^){n.i-^), (1) 



which describes fermions of spin one half which hop on a 
lattice with hopping amplitude Vij and interact on each 
site with a repulsion energy U. An interaction quench has 
so far been experimentally realized in the bosonic version 
of the Hubbard modelJ^ To describe the corresponding 
situation in the fermionic model, we allow for a time- 
dependent interaction U (t) in Eq. Q . 

Even after the mapping to a single-site model, the solu- 
tion of the Hubbard model within nonequilibrium DMFT 
requires the calculation of the time evolution of an inter- 
acting many-body system. In a previous publication^^ 
we employed a recently developed diagrammatic impu- 
rity solver— to compute the time evolution after an in- 
teraction quench over a wide parameter regime within 
DMFT. The numerical simulations confirmed an analyt- 
ical flow equation analysis for quenches to small U^^ 
which showed that in the limit U ^ the system is 
trapped in a nonthermal metastable intermediate state, 
a phenomenon known as prethermalization.— We identi- 
fied a similar trapping phenomenon for quenches to very 
large interactions. Most interestingly, these two prether- 
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malization regimes are separated by a well-defined "crit- 
ical" interaction U^ym where instead of a trapping in ei- 
ther of the nonthermal states a fast thermalization is ob- 
served. In Ref. [10 these qualitatively different regimes 
were demonstrated on the basis of an analysis of the mo- 
mentum distribution and the double occupancy. 

Relaxation to a thermal state is often difficult to es- 
tablish numerically because the time evolution must be 
studied on long timescales. In general thermalization 
is expected for interacting systems (in the weak sense 
that the expectation value of a large class of observables 
approaches the thermal expectation value in the long- 
time limit). In exactly solvable systems, however, ther- 
malization is often prevented by integrability j^^i^^'^^i^^ 
while its mechanism for nonintegrable systems is cur- 
rently under debatej ^'''^^i^^'^°i'^^'^^ The existence of two 
separate relaxation regimes is similar to what was found 
for Heisenberg spin chains^ and the one-dimensional 
Bose-Hubbard model.— 

The purpose of the present work is twofold: First, we 
want to explain in some detail the machinery behind our 
Quantum Monte Carlo (QMC) calculation of the Hub- 
bard model in nonequilibrium DMFT. We will briefly 
review the DMFT formalism (Sec. |TT| and the diagram- 
matic Monte Carlo method (Sec. IIIip . discuss some im- 
portant tricks which improve the efficiency of the Monte 
Carlo sampling, and then present in detail the solution 
of the DMFT self-consistency equations based on the ex- 
act equation of motion approach (Sec. IIV[) . The QMC 
solution of DMFT is finally used to discuss the validity 
of the nonequilibrium generalization of the iterated per- 
turbation theory (Sec. |V|. The second purpose of this 
paper is to further analyze the main finding of Ref. [20I . 
namely a fast thermalization after a quench from U ~ 
to Udyn, with additional data for the momentum distri- 
bution, the spectral function, and the optical conduc- 
tivity (Sec. IVip . In particular we find that at J7dyn the 
retarded non-equilibrium Green function relaxes to the 
appropriate equilibrium Green function within the nu- 
merical accuracy, establishing thermalization of all one- 
particle quantities that can be calculated from it. 



II. NONEQUILIBRIUM DMFT 

A. Contour-ordered Green functions 

In the following section we set up the framework for 
the investigation of a rather general class of nonequilib- 
rium situations. We assume that the system of interest 
is initially prepared in thermal equilibrium. For times 
i > it is then acted on by some perturbation, but 
there is no coupling to external heat or particle reservoirs. 
Technically, this setup implies that the time evolution 
is unitary and captured by a time-dependent Hamilto- 
nian, but all results must be averaged over initial states 
according to the grand-canonical density matrix po 
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FIG. 1: The L-shaped contour C for the description of tran- 
sient nonequilibrium states with initial state density matrix 
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PH(o)^ Arrows show a possible Monte Carlo configura- 
tion corresponding to perturbation order n — 7 and n+ = 3, 
n_ = 2, = 2 (cf. Sec. HIT)) . 



ventional approach to this kind of nonequilibrium situ- 
ation within many-body theory is based on the use of 
contour-ordered Keldysh Green functions,— 



(2) 



where the time arguments t and t' lie on the L-shaped 
contour C that runs from to some time tmax (i-e-, the 
largest time of interest) on the real time axis, back to 0, 
and finally to — i/? along the imaginary time axis (Fig.[T]). 
Here and in the following, operators with hat are in 
Heisenberg notation with respect to the time-dependent 
Hamiltonian [on the imaginary branch, H = H{0), and 
c(-ir) = e^^(°)ce-^^(°)], and (•) = Tr[po-] is the expec- 
tation value taken in the initial equilibrium state. The 
contour-ordering Tc exchanges the order of two operators 
A{ti) and B{t2) in a product A{ti)B{t2) if and only if t2 
appears later on the contour than ti, with an additional 
minus sign if the exchange involves an odd number of 
Fermi operators. The order of time arguments along C is 
indicated by the arrow in Fig. [TI which points from "ear- 
lier" to "later" times. Contour-ordered Green functions 
were first introduced by Keldysh^ in order to generalize 
Wick's theorem and diagrammatic perturbation theory 
to nonequilibrium physics. The extension of the origi- 
nal Keldysh formalism to the L-shaped contour C, which 
has numerous applications in nonequilibrium many-body 
theory,— becomes important whenever correlations be- 
tween the initial state at i = and time t > cannot be 
neglected 

The contour-ordered Green function ^ is related to 
a number of real and imaginary-time Green functions, 
which we list in the following paragraph for later refer- 
ence. When both time arguments are on the imaginary 
branch, Eq. ^ reduces to the Matsubara Green function 
of the initial equilibrium state, 



Gaa' ('^) '^') = Gaa' (-^T, -Zt'). 



(3) 



/Tr[e-''3^(o)] at temperature T = The con- 



Because the Hamiltonian is constant on the vertical 
branch of C and commutes with the initial state density 
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matrix, G"^, (r, r') is translationally invariant in imagi- 
nary time, such that we can introduce the usual Matsub- 
ara frequency representation, 



(4a) 



/3 



(?"(iw„) = -iJ dTe*'^"^G"(T,0). (4b) 



On the other hand, when both time arguments are real, 
one obtains the lesser, retarded, and advanced Green 
functions, 

G<^,{t,t') = G^^,{t+,t'_)^i{cl{t')c^{t)), (5) 
GL,(<,<') = e(t-t')[Gaa'(i-,t'+)-G„„Ki+,i-)] 

= -zeit-mcUt'),Ut)}), (6) 



Gto.,{t,t') = e(t'-0[Gaa'(i+,t'_)-G„a. (*-,*'+)] 



= ^e{t'-t){{cl{t'),c^{t)}). 



(7) 



The subscript of each real time-argument indicates 
whether it is on the upper {+) or lower (— ) real-time 
branch of C. The lesser Green function is related to the 
occupation of states a, to which its imaginary part re- 
duces for t — t' and a = a' . On the other hand, the 
retarded and advanced Green function arc related to the 
spectral function, which will be discussed in more detail 
in Sec. IVIl In addition to the real and imaginary time 
Green functions, the Green functions 



^aa'i^T''') = GQ,Q'(t±, — it), 
G|^a'(T, <) = Gaa'i-iT,t±). 



(8a) 
(8b) 



with mixed time arguments encode the correlations be- 
tween the initial state and times t > 0. 

It follows from the cyclic property of the trace and the 
definition of the contour-ordering that the Green function 
(PI satisfies an antiperiodic boundary condition on C in 
both time-arguments, 

Gaa'(0+,O = -Gaa'(-i/3,<'), (9a) 
G„a.(i,0+) = -G„„.(i,-i/?). (9b) 

This boundary condition holds for all contour functions 
in this text, including those which have no simple defini- 
tion in terms of Heisenberg operators. Furthermore, the 
Green function ([2|) satisfies the hcrmitian symmetry 

G^^,{t,t') = Gt,At',tr (10a) 

GL'M ^ -G<,Jt\tr (10b) 

G:„.(t,r) = G:,„(/3-r,0*, (10c) 

which will be used frequently in the following. 

B. Dynamical mean-field theory 

In equilibrium DMFT local correlation functions are 
obtained from a single-site impurity model subject to 



a self-consistency conditionii The mapping of the lat- 
tice problem ^ onto the single-site problem is formally 
achieved by integrating out all lattice sites apart from 
one. A straightforward reformulation of this mapping 
for Green functions on the Kcldysh contour— makes 
DMFT applicable to noncquilibrium problems. The 
single-site action is then given by 



5 = 5o + / dthioc{t) 



So 



= E 

<T=T,ic 



dtdt'cUt)A,{t,t')c,{t'), 



(11a) 
(lib) 



where J^dt — J^""''^dt^~J^""''^dt^^i Jq dr is the integral 
along C, 



/»ioc(i) = t/(<)(^T-i)K-5) 



(12) 



is the local interaction of the Hamiltonian, and Sq de- 
scribes the hybridization of the site with an environment 
that is determined sclf-consistcntly by the DMFT pro- 
cedure. In the following we consider only homogeneous 
paramagnetic phases, such that A(j does not depend on 
the lattice site or spin a. 

The local Green function for action (fTTj) is given by 



G„{t,t') = -i{c,{t)cUt'))s, 



(13) 



where operators without a hat arc in the interaction pic- 
ture with respect to fj.{nij + rt^), and the notation 



Tr[e-'3^("T+"i)Tc exp{-iS) ■ ■ ■ 
Tr[e-^'^("T+"i)Tc exp(-i5)] 



(14) 



is used. In general, the computation of Ga{t,t') is a 
complicated noncquilibrium many-body problem. For 
this reason, noncquilibrium DMFT has so far been ap- 
plied mostly to the Falicov-Kimball model, where the 
single-site problem can be reduced to a quadratic one and 
thus becomes exactly solvable either numericallyi2*i^ or 
analyticallyji^ In the present paper, just as in Ref. [20l . 
we investigate the Hubbard model and solve the single- 
site problem using the weak-coupling continuous time 
Monte Carlo algorithm^ which will be described below 
(SecHni. 

The local self-energy is then defined by the Dyson 
equation 



[{Go,l-^a)*G,]{t,t')^Sc{t,t'), 



(15) 



where the noninteracting ([/ = 0) single-site Green func- 
tion and its inverse are given by 

Go,a(t,t') = ~i{c,it)4{t'))s„, (16a) 
Gq-i (t,t') = Sc{t,t'){^^t+^i) - Aa{t,t'), (16b) 

respectively. Here we introduced the notation [a* 6](t, t') 
= dt a{t,t) b{t,t') for the convolution of two contour 
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functions, and the contour-delta function 5c(t,t') is de- 
fined such that 



dtf{t)Sc{t,t)=f{t) 



(17) 



for any contour function /, i.e., 5c{t,t') ~ ±S{t — t') if 
t and t' both on the upper or lower real branch of C, 
and (5c(— ir, —it') = id{T — t') for time arguments on the 
vertical branch. Both the Dyson equation (jisp and the 
corresponding equation 



(18) 



for Gq.ct are inhomogeneous integro-differential equations 
on the contour C when Eq. p6ap is inserted. They have a 
unique solution because Go,cr and G satisfy the boundary 
condition ^ . The solution of such integral equations on 
C is discussed in detail in Sec. HVl 

In order to determine the hybridization function 
Acr{t, t') one must equate the self-energy Ti^it, t') and the 
Green function Ga{t,t') of the single-site problem with 
the local self-energy Yijj„{t,t') and the local Green func- 
tion Gjja-{t,t') of the lattice problem at the given site j, 
respectively, 

Gj,^it,t') = GAt,t'), E,;,<,(i,i') = 5„E,(i,i')- (19) 

The latter two are related by the lattice Dyson equation, 

(^a^ + ek)Gka {t, t') - [S, * Gk.] (t, t') = 5c {t, t'), (20) 

which is stated here for the homogeneous case after 
Fourier transform with respect to lattice sites. In 
Eq. jSOl), 



Gfc,(t,t') = -*(TcCfc.(t)cL(i')) 



(21) 



is the momentum-resolved lattice Green function. (For 
a Bravais lattice, k are quasimomenta and are band 
energies, but more generally, and efc are eigenvectors 
and eigenvalues of the hopping matrix Vij, respectively.) 
The local Green function is given by the momentum sum 



(22) 



which closes the self-consistency. 

In the present paper we consider the case of a time- 
dependent interaction but no external fields. The hop- 
ping matrix elements are then independent of time, and 
the fc-summation in Eq. (j22p can be reduced to an inte- 
gral over a single energy variable 



G,it,t')^ dep{e)G,^it,t'), 



(23) 



involving the Green function G^„{t,t') = Gkc!{t,t')\^^^, 
and the local density of states p{e) = \ 



Efc) at an arbitrary site j. For the case of a semielliptic 
density of states. 



2ttV 



(24) 



with quarter bandwidth V, which corresponds to nearest- 
neighbor hopping on the Bethe lattice2.i^ or a particular 
kind of long-range hopping on the hypercubie lattice,— 
one obtains a closed form expression for the Weiss fieldfi^ 



A„it,t') = V^G,{t,t'). 



(25) 



We will use this self-consistency equation for all results 
of this work, so that the solution of the DMFT equations 
is achieved by iteration of Eqs. and ([^5]) . 



III. REAL-TIME MONTE CARLO METHOD 

The real-time evolution of the impurity model can be 
computed using the weak-coupling diagrammatic Monte 
Carlo method. More specifically, we employed the real- 
time version of the continuous-time auxiliary field al- 
gorithm (CTAUX)^ which was discussed in detail in 
Ref. [2l|. In the following section we present the im- 
plementation of this algorithm on the L-shaped contour 
(Fig. [1]) and then discuss some technical aspects which 
improve the efficiency of the method to the point where 
relevant timescales even in the strong-coupling regime 
become accessible. 

We start by expressing the partition function of the 
initial state as 

Z = Tr[e-'^'"("T+"i)7;,e-"5] 

^ g-/3(7/4-i/£,dtfe(t)rp^|-g-/3M(«i-+«i)2i^g-t5fcj (26) 

with Sk=S-J^ dt{k{t) + U/A) and k{t) ^ 0. This (pos- 
sibly time-dependent) shift in the action is introduced so 
that the interaction term can be decoupled using auxil- 
iary Ising spin variables according to^ 



hioc{t) + k{t) + [//4 



cosh(7(t)) 



k{t)-U{n^n^^{n^+ny)/2)) 
k{t)/2 e'^(*)"("T-"i)^ 



s=-l4 
U/{2k{t)). 



(27) 



Expansion of e^''^'^ in powers of (/iioc(0 ^ k{t) — U/4) 
and subsequent auxiliary field decomposition leads to 
an expression of the partition function as a sum over 
all possible Ising spin configurations on the contour 
C. The weight of the Monte Carlo configurations 
{(ti, si), (^2, S2), . . . (t„, s„)} (see illustration in Fig. [T|) 
is obtained by evaluating the trace of the remaining non- 
interacting problem. 
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si), . . . (i„, s„)}) - {ik{h)dt/2) . . . {ik{t,,^)dt/2){^ik{tn^+,)dt/2) . . . {-ik{t,,^+^_)dt/2) 

x{k{tn^+n_+l)dT /2) . . . {k{tn^+n_+np )dTl2)\{dctN-\ (28) 



-1 



(29) 



r 



Here n± and denotes the number of Ising spins on the 
three branches ofC, Gq^o- is the (n++n_+n_g)x (n++n_ + 
n/j) matrix of bath Green functions (|16ap evaluated at the 
time arguments defined by the Ising spins, and e^" ^ 
diag(e''(*i)"i'", . . . , eT^*")""""). A Monte Carlo sampling of 
all possible spin configurations can then be implemented 
based on the absolute value of these weights. 

The contribution of a specific configuration c to the 
Green function is given by 

G^(t,i') = Go,.(t,t') 

n 

+ I Go,.(t,i«)[(e^- - l)A^,],,,Go,.(i,,t'), 

(30) 

so the Green function G is obtained as the Monte Carlo 
average of G'^. 

The sign problem in this method grows exponentially 
with the average perturbation order on the real-time por- 
tion of C To reach long times or strong interactions, it 
is therefore important to reduce this perturbation order 
as much as possible. In the particle-hole symmetric case, 
i.e., at half-filling and for a symmetric density of states, 
the parameter k{t) of the algorithm can be chosen such 
that only even perturbation orders appear in the expan- 
sion. In fact, for 



we have j{t) = in, e^*-*-**"^ = —1 and hence the spin 
degree of freedom effectively disappears. The algorithm 
then becomes the real-time version of Rubtsov's weak- 
coupling method^ for the particle-hole symmetric inter- 
action term h\oc{t) = U{n^ — ^){ni — ^) with weight 



,(ti,...,t„) = {~iUdt)"+ (iUdt)"- i-Udr)'''' 



Y[dct (iG, 



'O.a - -^I 



(32) 



k{t) = -U/A 



(31) 



(For a detailed discussion of the equivalence between the 
Rubtsov and CTAUX methods for the Anderson impurity 
model, see Ref . fisl) . The above choice of k{t) requires the 
implementation of Monte Carlo updates which change 
the perturbation order from n to n±2. We found, how- 
ever, that the odd perturbation orders are continuously 
suppressed as k{t) approaches —U/A, so one may as well 
choose k{t) = —U/A + S (with small S) in combination 
with rank one updates. 

The efficiency of the Green function measurement can 
be improved dramatically by the following simple tricks. 
First, we rewrite Eq. ([50]) as 



G^(t,t') = Go,.(t,t')+ [dsi [ds2GoAt.si)0 Yl Hsi,U)[ie^' - l)N,],^,5cis2,t,)) Go,.(s2,0, (33) 



r 



where the variables si and S2 run over the contour C, and cost. Furthermore, comparison of Eq. ([55)1 to the Dyson 
{■)mc denotes the Monte Carlo averaging. It is therefore equation shows that X is related to the self-energy 



sufficient to accumulate the impurity system T-matrix 



Xa{si,S2) = 51 ^c{si,k)[{e^'-'^)N„],^jScis2,t,) 



by 



X * Go = S * G, 



(35) 



so the measurement of X allows to extract E as explained 



(34) in Section [IV] C 



as mentioned in Ref. k2- While the measurement of X on 



Further improvements are possible. Assuming that the 



some fine grid introduces discretization errors, these can perturbation order on the real-time branch is non-zero, it 
be made negligibly small at essentially no computational follows from Eq. ([25[l that the weight of the Monte Carlo 
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configuration changes sign if the last spin (corresponding 
to the largest time argument) is shifted from the forward 
contour to the backward contour or vice versa. Since the 
absolute value of the weight does not change, these two 
configurations will be generated with equal probability. 
As a result, all terms in Eq. (j34|) which do not involve the 
last operator on the contour will cancel on average. It is 
therefore more efficient to accumulate only the contribu- 
tions to Eq. (p4|) from those pairs in which either 
i or j corresponds to the last operator on the real-time 
branch. (If all spins sit on the imaginary-time branch, 
no such simplification is possible.) We also note that the 
error bars on measurements can be substantially reduced 
by appropriate symmetrizations of the real and imagi- 
nary parts of X (symmetry lines si = <max, S2 = imax, 

Sl = S2)- 

IV. WEAK-COUPLING CTQMC + DMFT 

To use the weak-coupling CTQMC as an impurity 
solver within DMFT, we iterate the following two steps 
until convergence: (i) The local Green function Ga{t,t') 
is determined in CTQMC [Eq. (p3|) ]. using the noninter- 
acting bath Green function G'o,cr(i,i') as input, and (ii), 
Go,cr{t,t') is determined from its inverse (jl6b[) . using the 
QMC output Gait.t') and the self-consistency Eq. (^5)) . 
We start the iteration from an initial guess for Gn^cr{t, <'), 
for which we usually take the noninteracting equilibrium 
Green function, 

Gl^{t,t') = 1 1 dep(e)e-(*'-*)[/(e) - Qc{t,t% (36) 

where &c{t, t') = 1 if t is later on the contour than t' and 
otherwise zero. 

In this section we describe in detail how Go.cr is deter- 
mined from Act (Sec. IIV A[) . how the self-energy is cal- 
culated from the impurity correlation function after 
convergence of the DMFT iteration (Sec. IIV B|) . and how 
one finally obtains expectation values of various observ- 
ables of the lattice system (Sec. IIV C[) . Furthermore, we 
introduce a real frequency representation which is needed 
to efficiently treat the case of zero temperature on the L- 
shaped contour (Sec. lIVPp . and combine this with the 
weak-coupling impurity solver for the case of a noninter- 
acting initial state. 

A. Integral equations on the contour C 

Within nonequilibrium DMFT one must frequently 
solve equations on C of the type 

[idt - hit)]Yit,t') - [K * Y]it,t') = Scit,t') (37) 

with a known integral kernel K(t,t'). The solution 
Y{t,t') is unique when the antiperiodic boundary con- 
dition ^ is imposed on Y{t,t'). For example, both 



Eq. (jlSp for the noninteracting bath Green function 
Go.ij{t, t') and Eq. ([20]) for the momentum-resolved Green 
function have this form. 

By choosing a suitable discretization of the contour C, 
Eq. (|37p can in principle be reduced to the inversion of a 
matrix whose dimension is given by the number of mesh 
points along In the following we pursue a different 
approach, where both Y and K in Eq. (|37|) are first rep- 
resented in terms of their respective real and imaginary 
time components ([H])-®, and separate integral equations 
(which are similar to the Kadanoff-Baym equations^) 
are solved for each component. Although this proce- 
dure may seem rather cumbersome compared to direct 
contour discretization, it has several advantages: (i) It 
is straightforward to incorporate the hermitian symme- 
try pUj) which is satisfied by both the local self-energy 
and the hybridization function Aa-. (ii) The resulting 
equations are Volterra type integro-differential equations, 
for which highly stable and accurate algorithms can be 
found in the literature;^ and which remain causal even 
when they are approximated numerically. Finally, (iii), 
the real-frequency representation which we introduce in 
Sec. IIVDI to handle initial states at zero temperature is 
based on this approach. 

In the following we assume that Y and K satisfy the 
hermitian symmetry (|10p , such that it is sufficient to de- 
termine the Matsubara, retarded, mixed , and lesser 
components of Y. Corresponding components of the con- 
volution K*Y in Eq. ((37|) are obtained from the Langreth 
rules which follow directly from the definitions Q-® 
and the definition of the contour integral. By taking the 
Matsubara component ^ of Eq. ([57)1 we obtain 

i-dr-h)Y''iT,T')+i J dfK"{T,f)Y"{f,T') = 


i(5(T-r'), (38) 

where h = ft,(0) is constant on the imaginary branch. 
This equation must be augmented with an antiperiodic 
boundary condition Y"[Q,t') = -Y^'[f3,T') which fol- 
lows from Eq. When we assume that the kernel if" 
has the Matsubara frequency representation ([4]), it fol- 
lows that the solution F"(r, r') is of the same form, with 

y"{iu;n) = [iuJn-h-k"{iu)Y\ (39) 

As required by causality, Y" thus turns out to depend 
only on the initial equilibrium state, independent of the 
subsequent perturbation of the system. 

In a similar fashion, the retarded component ^ of 
Eq. ([371) is given by 

t 

[idt - h{t)\Y^{t,t') - j dtK^{t,t)Y^{t,t') = 5{t - t'). 

(40) 
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Because Y^{t,t') vanishes for t < t' hy definition 
[cf . Eq. ^] , integration over the (5-function yields 



(41) 



One can thus restrict the sohition of Eq. ([1D|) to t > t', 
drop the (5-function on the right-hand side and instead 
impose (|4ip as an initial condition. 

The limits of the integral in (|40|) take into account 
that retarded functions vanish for t > t' . This fact turns 
Eq. (|40|) into a Voltcrra equation of second kind,— i.e., 
the derivative at time t is determined by the kernel and 
the function at earlier times only. The numerical solution 
of this type of equations is analogous to the solution of 
ordinary differential equations4^ 

For the Green functions ([5]) with mixed time argu- 
ments, Eq. ([57]) reads 



[idt - h[t)]Y^[t,T) - j dtK^{t,t)Y^{t,T) 







^-i J dfK^it, f)y " (f , r) . (42) 



We assume that Y is continuous on C (which is true 
if neither K{t,t') nor h{t) are singular at t = 0), such 
that Eq. (|^^ must be solved with the initial condition 
y^(0, r) = y"(0, r). For given r, Eq. (gH) is an inhomo- 
geneous Volterra integro-differential equation, for which 
only known functions [cf. Eq. (|39p ] enter the source term 
on the right-hand side. 

A third and last Volterra integral equation can be de- 
rived for the lesser component (O, 



[idt - h{t)]Y<{t,t') - J dtK^{t,t)Y<{t,t') = 



P t' 
-iJdfK^{t,f)Y"{f,t')+ jdtK<{t,t)Y^{t,t'). (43) 



Due to the symmetry (fTO|) it is sufficient to solve this 
equation for t < t' , with the initial condition F^(0, t') — 
-y^(/3, t'). The latter follows from Eq. © and the con- 
tinuity of Y along C. The functions Y^ and Y^ which 
enter the source term of Eq. (|43|) on the right-hand side 
can be obtained from the previous solution of Eqs. (PD|) 
and (im, and the symmetry pI7|) . The successive so- 
lution of Eqs. ((all), (go]), dMl), and ([43]) completes the 
determination of the contour function Y . 



B. Determination of the self-energy 

The impurity self-energy can be obtained from the cor- 
relation function via Eq. (j35p . By comparison of 



the Dyson equation ((T5|) in integral form, Ga = Go,a 
+ G„ * Go,,,, with Eq. i.e., G„ = Go.a + 

Go,a * Xa * Go,cr, wc find the relation 



(44) 



This equation is very similar to Eq. (j37p . with unknown 
y = S, kernel K = X^r * Gq.o-, h ~ 1, and without the 
differential term. The solution of (jH]) is thus analogous 
to Eq. (|37p . using a decomposition in terms of the com- 
ponents ([21)-©. The final equations read 



1 -I- A:"(itj„)' 



(45a) 



J:^it, t') + J dtK^it, t)j:^{t, t') = X^{t, t'), (45b) 
f 

t 

YT {t, t)+ J dtK^it, t)S" {t, t) = X^ {t, t)+ 



i j dfK^{t,f)T.'''{f,T) (45c) 



t 

^<it,t') + j dtK^{t,t)i:<{i,t') ^ x<{t,t')+ 



/3 t' 

i j dfK^ {t,f)Yr{f ,t')-jdt K<{t,T)T.\t,t'). (45d) 



Note that the kernel K = X^, * Go,cr does not satisfy the 
hermitian symmetry ([TOll . i.e., X^, * Go,o- ^ Go,^ * Xg-. 

We would like to remark that the self-energy can 
equally well be determined from the linear equation 



Xcr * GctO = So- * Go 



(46) 



However, Eqs. (gS)) are essentially Volterra integral equa- 
tions of the second kind, while Eq. ([46)1 leads to Volterra 
equations of the first kind, i.e., only the integral-term 
is present on the left-hand side. Because the numeri- 
cal solution of Volterra equations of the first kind tends 
to be unstable^ we prefer the solution of Eq. (144]) over 

Eq. dm). 



C. Expectation values of observables 

From the self-energy S one can directly compute the 
expectation values of observables of the lattice Hamilto- 
nian. In this section we let (• • •) denote the initial state 
expectation value at temperature T = 1//?, and operators 
with hat are in Heisenberg representation with respect to 
the Hubbard Hamiltonian ([1]) with time-dependent inter- 
action. The number of lattice sites will be denoted by L. 
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The particle number per site for spin a is given by the 
local Green function G„{t,t') 

nAt) ^jY^i^lmAt)) = -^G<{t,t), (47) 



provided that the state is homogeneous. Because n^{t) 
is conserved, the condition G^{t,t) — const provides a 
first test of the numerical accuracy. 

The occupation of the momentum states 

n{ek,t) EE (cLW4.W> = -*G<,(i,t), (48) 

is obtained from the momentum-resolved Green function 
Gka{t,t') = ^*(TcCfccr(t)c^^(t')). For a momentum in- 
dependent S^, n(ek,t) depends on momentum k only 
via the band-energy Cfc. The Green function Gka{t,t') is 
determined from the lattice Dyson equation (PO)) . whose 
solution is analogous to that of Eq. ([57]) . The kinetic 
energy per lattice site 



i?kin(t) = ;^^efc(cLWcfc(i)), 



(49a) 



is obtained from n{e, t) by replacing the fc-sum with an 
integral over the local density of states [Eq. ([24]) ]. 



£'kin(i) = / dep{e)en{e,t). 



(49b) 



Furthermore we are interested in the double occupa- 
tion per lattice site 



(50) 



and the interaction energy 

£;pot ^ U{t) J2 { («.T (i) - 5) ) (51) 

i 

= U{t)[d(t)-^{n^{t)+n^{t)) + l]. (52) 

To calculate this quantity we consider the equation of 
motion for the local lattice Green function Gjja, which 
reads 

[{G-^)ji * Gi,,]{t,t') ^ Sc{t,t') + U{t)V,At,t'), (53) 



{G-A),i{t,t') = 5c{t,t')[5fl{idt + ii) - t,i] 



(54) 



T,,{t,t') = -i{Tcha{t)[n,s{t) - \)cl{t'))- (55) 

Comparison with the lattice Dyson equation in real space 
yields 



U{t)Tj„{t,t') = [E,*G„,](t,t'), 



(56) 



because the self-energy is local and site- independent. 
Hence = Tj„ can be determined from quantities 
measured in the single-site problem [cf. Eq. ((46|) ]. and 
Eq. fSS]) implies 



for a homogeneous state. 

Finally we can compute the total energy from Eqs. 
and ([5T|). 

£^tot(0 =^ki„(<) + Spot(^)■ 



(58) 



This quantity must be constant when the Hamiltonian is 
time- independent, which provides a second test for the 
accuracy of the numerical solution. 



D. Real- frequency representation 

In this subsection we introduce a partial Fourier trans- 
form of the mixed components and , which will 
allow us to handle contour equations such as Eqs. (|37[) 
and (j44p in the limit of zero temperature without dealing 
explicitly with a contour of infinite length. 

We start from the Fourier series in the interval < 

T</?, 



(59a) 



in terms of fermionic Matsubara frequencies iuJn- This 
representation is now used within the solution of Eq. (j37p . 
In contrast to Eq. (|^^ . the corresponding equation for 
the mixed component 



{-dr - h)Y^ (r, t)+i df K"{t, f)Y^ (f , t) 



Z 

JdtK"{T,t)Y^{t,t) (60) 



dit)=-ir<{t,t) + ^nAt) 



(57) 



is not an initial value problem, but rather a boundary 
value problem on C: The boundary condition Y^ (/3, t) = 
-F"(0,t) - ^"^(0,*) follows from Eqs. Q and © and 
the continuity of the contour functions along C. Using 
the transformation ([59]) and Eq. ([39]) for the Matsubara 
component, Eq. (|60p becomes an explicit integral expres- 
sion for Y^ {iun, t), 

Y^{iuJn,t) = 

t 

J/"(iw„)[rA(0,t)+ JdtK^{iLUn,t)Y^{t,t)\. (61) 


In the following we assume that K'~{iuJn, t) can be con- 
tinued to complex frequencies z, such that K^{z,t) is 
analytic in the upper and lower complex half plane, re- 
spectively, and has a branch cut along the real frequency 
axis. For the Green function @ this property follows 
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from a Lehmann representation in terms of an eigenbasis 
{|n)}ofi?(0), 

r<^ I ■sr^{w,,+w„){n\ca\m){m\c^^,{t)\n) 



Z + En — En 



which has poles on the real axis only {w„ 



When Eq. ([M]) is continued to the real axis, we obtain 
two functions [lo'^ , t), where uj^ ~ ojHrj for 77 ^ 0+. 
In contrast to the equilibrium functions the two 

functions Y'~{uj^,t) are not simply related by complex 
conjugation. Matsubara summations are then trans- 
formed into integrals along the branch cut of [z, t) 
in the usual way. For example, the backtransformation 
(|59bp is given by 



y-(r,t) 



duj 



/(c.)e--[r^(c.-t)-r^(c.+,<)], (63) 



where /(w) = \/{e~^^ + 1) is the Fermi function. 

Using the real frequency representation for the mixed 
components, the first source term on the right hand side 
of Eq. (|43p can be rewritten in terms of the previously 
determined functions Y^ {uj^ ^t) 



/9 

JdfK-it,f)Y-(f,t) = j^^f{. 



uj] X 



[K^ {t, Lo+)Y- {uj+ ,t')-K^{t,LU-)Y- {u;~ ,<')]• (64) 

Here the Fourier transformation (j59ap with opposite sign 
for iujn is used in the second time argument for the ^ 
Green functions (I5al). such that 



(65) 



Y-{z,t) = -Y^it,z*r 



follows by symmetry (jlOcp . 

The above real-frequency representation can be used 
within DMFT whenever the impurity problem is solvable 
at zero temperature. This is the case for approximate 
analytical methods (Sec. |V|. It might also be of ad- 
vantage for arbitrary initial states in the Falicov-Kimball 
modelfi^ii^ where the solution of the impurity is based 
on the solution of equations of motion which have ex- 
actly the structure of Eq. ((57)1 . In the following section 
we present another application, namely nonequilibrium 
DMFT for the Hubbard model with a noninteracting ini- 
tial state. 



E. CTQMC and DMFT for noninteracting initial 
states 

The computational scheme for the solution of the 
nonequilibrium DMFT equations is represented in Fig. [2] 
for a semiclliptic density of states and a noninteract- 
ing initial state. Green functions Go-, Ga-o, and A^- sat- 
isfy the symmetry (|10p . such that they are represented 




FIG. 2: Computational scheme for the nonequilibrium 
DMFT using the self-consistency p5|) and noninteracting ini- 
tial states. Steps (i)-(iv) are explained in the text. 



by their Matsubara, retarded, and lesser component. 
The Matsubara Green functions are given by the equilib- 
rium (noninteracting) Green function 



de- 



(66) 
(67) 



where p{e) is given by Eq. ([24|) . and the last equahty 
holds due to the self-consistency (|25p . The mixed com- 



ponents , GJ7o , and are represented after the par- 
tial Fourier transform (j59p and analytical continuation 
by their value along the branch cut, i.e., for each Green 
fimction Y we keep two functions Y^ {uj^ ,t) on a fixed 
frequency mesh. 

The DMFT iteration is started from an initial guess 
In step (i) [cf. Fig. the Weiss field A{t,t') is 
computed from the closed self-consistency equation (|25p . 
This is then used to determine the noninteracting bath 
Green function Gao from its inverse (jl6bp . as explained 
in the previous subsection [step (ii) in Fig. [5] . The func- 
tion Gcro is the input for the calculation of the interacting 
bath Green function ([13]) using CTQMC [step (iii) in Fig. 
[2] . Because the initial state is noninteracting, the Monte 
Carlo simulation is restricted to the real-time branch 
of the contour, and only the real-time components G^ 
and G^ are obtained. However, the mixed component 
G'^{uj'^,t) can be reconstructed from these functions and 
the previous Weiss field A [step (iv)]: For this purpose 
consider the Dyson equation p^ . which has the form of 
Eq. ([37]) after the replacement if = A^, -|- Eg., F = G^, 
and h{t) = fjL. Hence G^(z,i) can be obtained from the 
integral (pTj) . making the same replacements. Because 
Ti„{t, t') is proportional to the interaction strengths U{t) 
and U(t'), we have S^(r, t) — Q for a noninteracting ini- 
tial state, with C/(-ir) = C/(0) = 0. Hence G^(w±,i) is 
given by 

G:(c.±,t) 

t 

= 5^(^±)[G^(0,t)+ JdtA-{uj^,t)G^it,t)], (68) 


where g"(w^) = ^nip{uj) [Eq. (|66p ]. Steps (i) through 
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FIG. 3: Potential energy Epot{t) [Eq. ([51) ]. kinetic energy 
Eun{t) [Eq. dMJ], and total energy Etot = Epot + -Bkin in the 
half-filled Hubbard model after an interaction quench from 
[/ = to [/ = 2 (a) and !7 = 5 (b). The initial state tem- 
perature is given by T = 0. The results were obtained with 
nonequilibrium DMFT for a semielliptic density of states (|24[) . 
using either CTQMC (symbols), IPT (dashed lines) or SPT 
(solid lines) to solve the single-site problem. The inset in the 
upper panel shows the second-oder diagram for E^. Lines 
represent Go,a for IPT [Eq. (glj], and G„ for SPT [Eq. JTO}]. 



(iv) are repeated until convergence, which is usually 
reached after not more that 15 iterations. 



V. COMPARISON TO ITERATED 
PERTURBATION THEORY 

In equilibrium DMFT, the so-called iterated perturba- 
tion theory (IPT)^i^ is frequently used as an approxi- 
mate but efficient method to solve the single-site prob- 
lem. Within IPT, the self-energy E^r is expanded up to 
second order in the interaction U. Although this is a 
weak-coupling expansion by construction, it is acciden- 
tally correct for the atomic limit of the half-filled Hub- 
bard model in equilibrium. In many aspects. IPT thus 



provides a reasonable interpolation between the two ex- 
act limits U = and ^ = 0. In particular, it qualita- 
tively reproduces the DMFT phase diagram and the Mott 
transition in the paramagnetic phase, although there are 
quantitative differences to numerically exact QMC re- 
sults. It is therefore interesting to see whether this ap- 
proximation performs similarly well when used to solve 
the single-site problem in nonequilibrium DMFT. 

In the following we restrict ourselves to the half-filled 
Hubbard model with time-dependent interaction U{t). 
The Hartree contribution to the self-energy (first-order 
diagram), which gives a shift of the chemical potential 
with respect to /i = 0, then vanishes, and the second- 
order contribution to the self-energy is given by a single 
diagram (inset in Fig. [3^), 



E|P*(t,i') = -Uit)Uit')Go,At,t')Go,ait',t)Go,a{t^t'). 



(69) 



This equation is easily incorporated into the DMFT self- 
consistency iteration by replacing step (iii) and (iv) in 
Fig. [5] with a solution of the Dyson equation (|15p for Go-, 
where Sg- is given by Eq. ((69l) . Equation ([15]) is solved 
numerically, as described in Sec. IIV Al 

In Fig.[3|we plot the potential energy Epot{t) [Eq. ([ST]) ]. 
the kinetic energy £'kin(i) [Eq- ([^51) ]- and total energy 
Etot = Epot + E]iin of the half-filled Hubbard model af- 
ter an interaction quench from the noninteracting ini- 
tial state at temperature T = 0. The hopping matrix 
elements correspond to a semielliptic density of states 
Eq. ([M]) with quarter bandwidth V = 1, and time is 
measured in units of h/V = 1. The numerically exact 
CTQMC results show a rapid relaxation of these quan- 
tities, which is discussed in detail below. As required 
by energy conservation, i?tot is constant within the nu- 
merical accuracy. IPT can reproduce these results rather 
accurately for small values of U (Fig. [3^). Already at in- 
termediate coupling, however, the results of CTQMC and 
IPT strongly deviate from each other (Fig. [3]d). In par- 
ticular, the total energy Etot is generally not conserved 
within IPT, such that the use of IPT as an approxima- 
tion for the intermediate- and strong-coupling regime be- 
comes highly questionable. In contrast to equilibrium 
DMFT, IPT does not provide a reasonable interpolation 
between weak- and strong-coupling regimes. 

This violation of energy conservation is cured by a sim- 
ple procedure. An expansion of up to finite order 
in terms of the noninteracting Green function is not a 
conserving approximation in the sense of Kadanoff and 
Baym42i^ However, the approximation becomes con- 
serving when Gq.ct in Eq. ([69)) is replaced by the full 
interacting Green function, 

i:Tit:t') = -U{t)U{t')G„{t,t')G^{t',t)Gs{t,t'). (70) 

The resulting self-consistent perturbation theory (SPT) 
is a truncation of the skeleton expansion for the self- 
energy, which can be derived from an approximation to 
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the Luttinger Ward-functional and is therefore conserv- 
ing. SPT is incorporated into the DMFT iteration by re- 
placing step (ii)-(iv) in Fig. [2] with a solution of the Dyson 
equation (fTS]) for Ga , where S^r is given by Eq. . Note 
that in this implementation is the SPT solution of the 
single-site problem for given A only after the DMFT it- 
eration is converged. 

When SPT is used instead of IPT as an approximate 
impurity solver, we find that E'tot is indeed constant with 
time (solid lines in Fig. [3|). However, SPT is not reliable 
at intermediate interaction strength either. For U = 5 
(Fig. [SJd), SPT predicts a monotonous relaxation of E'pot 
and i?kin, while the numerically exact QMC yields oscil- 
lations which are an important feature of the dynamics in 
the Hubbard model at strong coupling. For weak inter- 
actions, SPT performs slightly better, but in this param- 
eter regime it is worse that the IPT solution (Fig. [3^). 
The fact that IPT approximates the exact numerical solu- 
tion better than SPT is already known from equilibrium 
DMFT. 



VI. RESULTS 

In the remainder of this paper we present additional 
numerical results for the interaction quench in the Hub- 
bard model in nonequilibrium DMFT, building on our 
previous work (Ref. 120 ). The system is assumed to be in 
the noninteracting ground state before time t = 0, when 
the interaction is abruptly switched to a positive value 
U. We consider only homogeneous nonmagnetic states 
at half-filling (n^ = = i). Hopping matrix elements 
arc chosen such that the density of states is of semiellip- 
tic shape Eq. ((24)l . and the quarter bandwidth y = 1 is 
set as energy unit, so that time is measured in units of 
H/V=l. 

The time evolution of various thermodynamic quanti- 
ties after this interaction quench was already discussed in 
Ref. [20I . After some preliminary remarks on the effective 
temperature after a quench fSec lVl'S]) we will briefly re- 
state the basic conclusions of the latter publication and 
substantiate them with additional data (Sec. IVIB|) . We 
then turn to a characterization of the relaxing state in 
terms of dynamical quantities, i.e., the spectral function 
fSec lVICj) . and the optical conductivity (Sec. IVI P]) . 

A. Excitation after an interaction quench 

An important information on the state of the system 
after the interaction quench is its excitation energy with 
respect to the ground state. Because the system is as- 
sumed to be isolated from the environment, the total en- 
ergy is conserved after the quench and its value follows 
from the expectation values of the Hamiltonian in the 
initial state immediately before the quench. The energy 
corresponds to an effective temperature Toff, i-e., the tem- 
perature of the unique thermal equilibrium state which 




FIG. 4: Momentum distribution n(e, t) after an interaction 
quench in the Hubbard model from the noninteracting ground 
state to interaction [/ = 2 (a), [/ = 3.3 (b), and U — 5 (c). 

has the same total energy [Eq. 

Tr He^^^^'=" 
EMt) = i?tot(0+) = ^^^^H/n^^ ■ (71) 

(An analogously deflned effective chemical potential is 
fixed to /ioff = by particlc-holc symmetry.) For the 
quench in the Hubbard model we compute T^s by a nu- 
merical solution of Eq. (j71|) . Thermal equilibrium expec- 
tation values of static quantities are obtained from equi- 
librium DMFT, using QMC as impurity solver. For the 
quenches discussed below, Teg is of the same order as the 
hopping strength, which is far above the Mott transition 
endpoint in thermal equilibrium. 

If the system reaches a thermal equilibrium state a 
sufficiently long time after the quench, the temperature 
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of this state is given by Toff. Below we thus compare 
expectation values of observables after the quench with 
thermal equilibrium expectation values at T = Tcff. All 
static quantities in thermal equilibrium are directly com- 
puted within equilibrium DMFT. The computation of 
dynamical quantities such as the spectral function and 
the optical conductivity, however, would require an an- 
alytical continuation from Matsubara frequencies to real 
frequencies, which is not accurate enough at large fre- 
quencies and high temperature to allow for a quantitative 
comparison. We therefore use nonequilibrium DMFT to 
obtain real-time Green functions and the real-time op- 
tical conductivity in thermal equilibrium directly in the 
time domain. This calculation is equivalent to an "inter- 
action quench" in which the value of U is not changed 
and the initial state is at finite temperature T = Tcff- In 
contrast to the quench from the nonintcracting state, it 
is done on the L-shapcd contour, and wc do not use the 
tricks which are discussed in Sec. IIVEI The maximum 
times that are accessible in this way are comparable to 
the times which are accessible an the interaction quench 
from the noninteracting initial state. 



B. Relaxation after an interaction quench 

The time evolution after an interaction quench in 
the Hubbard model depends on the parameter U in a 
very sensitive manner. To illustrate the qualitatively 
different relaxation behavior in the weak, strong, and 
intermediate-coupling regime wc plot the momentum dis- 
tribution n(e,t) [Eq. gH)] for three values of U (Fig.H]). 
In all three cases the magnitude of the discontinuity 
An{t) = hm^^o+ ["-("'7: ~ '^(^)^)] a-t the Fermi energy 
decreases with time. Note that An{t) remains finite for 
a finite time after the quench; for the present case of a 
local self-energy this is due to the fact that An{t) is di- 
rectly related to the retarded Green function at e = 0^ 
Because a discontinuity in the momentum distribution of 
a Fermi liquid in thermal equilibrium can exist only at 
zero temperature, while on the other hand, a quenched 
system is always excited with respect to the ground state, 
the existence of a finite jump An{t) clearly indicates that 
the system is not yet fully thermalized. The size of the 
discontinuity is thus well suited to characterize the relax- 
ation after the quench. 

In the weak-coupling regime (Fig. [4^), n{e,t) rapidly 
evolves towards a distribution (t < 2 in Fig. 3^) which 
is not yet thermalized, but changes only slowly in time. 
This emergence of long-lived nonthermal states is an ex- 
ample of prethermalization^^ which is observed in a wide 
range of classical and quantum systems.— As shown by 
Moeckel and Kehrein^^ the nonthermal state remains 
stable for all times within second order unitary perturba- 
tion theory in U /V, i.e., higher-order corrections become 
effective only on the long timescale /U^. In the limit 
of infinite dimensions their weak-coupling result for the 
transient behavior towards the prethcrmalization plateau 
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FIG. 5: Comparison of the momentum distribution n{e, t) 
for fixed time t after the quench (symbols) to the momentum 
distribution in thermal equilibrium at the effective temper- 
ature Toft [cf. Eq. (|71|l ] (solid lines). Interaction parameters 
are [/ = 2 (a), U = 3.3 (b), and [/ = 5 (c). 



has the form 



«pcrt(e,t) =n(e)-4f/2T(e,t), 

°° sm^{E-e)t/2 



F{e,t) 



dE- 



(72) 
(73) 



— OO 

J^{E) = Jde[ Jde'^ Jdei S{e[ + e'^ - f-i - E) x 

P(ei)p(4)p(ei) [n{£)niei){l - n{e[)){l - nie'^)) 
-(l-n(e))(l-n(ei))n(e;)n(ei)]. (74) 



For a half-filled band and a symmetric density of states, 
p{e) = p(— e), wc obtain 



T(e,i) 



sgn(e) 



Z 

Jds{t~ s)Rc[R{sfe"'^'^], (75) 



where R{s) = /de e(-e) p(e) e*". This yields An(t) 
and also d{t) by using the energy conservation after the 
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oscillations is not fully accessible within CTQMC due 
to the dynamical sign problem. However, our results 
show that n{e,t) oscillates around a nonthermal distri- 
bution (Fig. [S}:). This behavior, which is analogous to 
prcthcrmalization at weak-coupling, is similar to what 
was found for the double occupation i.e., a decay 

on the timescale 1 /V to oscillations around a nonthermal 
value which does not change on much longer timcscalcs. 

The interaction quench to ?7 = 3.31/ is characterized 
by a rapid thermalization of the momentum distribution 
(Figs. HJd andlSb), without signatures of cither collapse 
and revival oscillations or a prethermalization plateau 
in n{e,t). Numerically we cannot detect a finite width 
to the crossover regime between the weak- and strong- 
coupling behavior, which indicates that there is a single 
point U — Udyn ~ 3.2V which marks a dynamical transi- 
tion in the Hubbard model.— A further investigation of 
this phenomenon and its relation to the Mott transition 
in equilibrium will require a systematic analysis of inter- 
action quenches which start from a wide range of initial 
states other than the noninteracting ground state. This 
is left to a future publication. In the following wc turn 
to a different question and investigate to what extent the 
rapid thermalization close to U = C/dyn, the oscillations 
at J7 > Udyn, and the prcthcrmalization at C/ < C/dyn 
show up in various dynamical quantities of the Hubbard 
model. 



C. Spectral function 



FIG. 6: Approach of the prethermalized state at weak- 
coupling and subsequent relaxation towards the thermal state, 
(a) Discontinuity at the Fermi surface, (b) Double occupa- 
tion. Solid lines: weak-coupling results [Eq. (|76|) -(I77 | )]. 

quench. 

t 

Anpert(0 = l-4U^ Jds{t~ s) Re[R{sf] , (76) 



t 

dpcrt (t) = ^ - 2C/ Jds Im [Risf]. (77) 



Numerical evaluations of these functions are plotted and 
compared to our DMFT results in Fig. [H] for the semicl- 
liptic density of states ([24|l with V = 1. Regarding the 
transient behavior and the prcthcrmalization plateau wc 
find very good agreement for ?7 < 1. Interestingly the 
prcthcrmalization plateau of An{t) is almost correctly 
predicted by the weak-coupling results even for [/ < 2. 
For larger times the system relaxes further towards the 
thermal value. 

In the strong-coupling regime (Fig. |4h), the relaxation 
is dominated by damped collapse and revival oscillations 
of approximate periodicity 2t:/U. The decay of these 



Important information about a correlated system out 
of equilibrium cannot only be obtained from thermody- 
namic quantities, but also from the dynamical response 
of the system to certain external perturbations, which 
can be computed from various real-time correlation func- 
tions. In the following subsection we discuss the time 
evolution of the local Green function Ga{t,t') = G{t,t') 
in the paramagnetic phase of the Hubbard model after 
a quench from the noninteracting ground state to finite 
interaction U. For this purpose we introduce the partial 
Fourier transform 

G^'<{uj,t) = Jdse"^'G^'<{t + s,t) (78) 

of the retarded and lesser Green function, and the spec- 
tral function A{u;,t) = -(I/tt) ImG^(w iO, i). The 
spectrum turns out to be a useful representation of the 
nonequilibrium Green function, although it lacks a direct 
relation to the "distribution function" G^ (w, t) and thus 
docs not have the same significance as in the equilibrium 
case. [In equilibrium one has G'^{uj) = 2TriA{uj)f{uj).] 

Before discussing the results we have to mention a tech- 
nicality, which arises from the restriction of the Monte 
Carlo simulations to relatively small times t < imax- In 
practice, the integration range in Eq. ([78|) must be cut off 
at Sniax = imax ~ t, leading to artificial oscillations at fre- 
quency 1/sinax- To reduce this effect in a controlled way 
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FIG. 7: (a) Local Green function G^{t + s,t) for an interac- 
tion quench in the Hubbard model to U = 3.3 (slightly above 
f^dyn). The function is purely imaginary due to particle-hole 
symmetry. The solid black line is the Green function in the 
thermal equihbrium state {U = 3.3, TeH = 1/13 = 0.84). The 
dotted line {U = 0) is the retarded Green function in the 
noninteracting initial state, (b) Spectral function A{uj,t) = 
— (l/7r)ImG^(a;,t) for the same parameters as in the upper 
panel. The dotted line and the line labelled T^s = 0.84 are 
the semielliptic density of states (|24|) of the initial state and 
the thermal equilibrium spectrum at temperature T^s = 0.84, 
respectively. Spectra are obtained from Fourier transforma- 
tion of real-time quantities, and the Fourier integral (|78|l is 
cut off at Smax = 3.5 with an additional Gaussian factor (see 
text). The corresponding kernel [Eq. (|79p . k, = 0.1] is shown 
as thin solid line. 



we introduce an additional Gaussian factor exp(— s^k) in 
the integral ([75]) . The resulting expression amounts to 
a convolution of the true Fourier transform (fmax = oo) 
with the kernel 



K, Sjnax) 



Smax 

— J ds cxp{iujs — s'^k). 



(79) 



A suitable choice of the parameter k can in some cases 
suppress the oscillations without washing out important 



spectral features, and a comparison with a known equi- 
librium spectrum is always possible without loss of in- 
formation after convolution of the latter with the same 
kernel. 

In Fig. [7] we plot G^ {t + s, t) and A{io, t) for a quench 
to interaction [/ = 3. The spectrum A{ijj,t) differs from 
the initial semielliptic density of states for all times t > 0, 
because the choice of the Fourier transform in Eq. (fTB]) 
implies that the initial equilibrium Green function does 
not enter the definition of A{u!,t) for t > 0. Note that 
this would be different for the common definition of the 
Fourier transform at constant average time {t + t')/2^ 
Within numerical accuracy, both G^ {t + s, t) and A{uj, t) 
become time- (<)-independent for t > \/V. This timescale 
is comparable to the relaxation time of the double oc- 
cupation and the momentum distribution at /7 ~ 3.3 
(Fig.Sb). 

An important interpretation of the finite relaxation 
time in A{ijj, t) can be inferred directly from the definition 
of the Green function. According to Eq. ([6]), G^ {t + s, t) 
is related to the survival amplitude of local single-particle 
excitations which arc created at time t and destroyed at 
later time t -\- s. The decay of such an excitation de- 
pends on both the Hamiltonian, which defines the pos- 
sible scattering mechanisms, and the quantum state of 
those particles which act as scatterers. While the Hamil- 
tonian changes abruptly at t = 0, the latter evolves with 
time, leading to the finite relaxation time of A(ijj^t). In 
contrast, A{uj,t) would be constant immediately after a 
quench in a noninteracting system, because the anticom- 
mutator in Eq. ([S]) is a c-number for a quadratic Hamil- 
tonian. Wc can thus conclude that the finite relaxation 
time observed in Fig. [7] is a true many-body effect, in 
analogy to the well-known fact that equilibrium spectra 
depend on temperature only for interacting systems. 

To characterize the final state after the relaxation, its 
spectrum should be compared to the equilibrium spec- 
trum of a correlated metal at rather high temperature. 
In fact, A{ijj,t) is strongly modified with respect to the 
semielliptic density of states, with precursors of the Hub- 
bard bands around lo = ±2. The fact that the spectrum 
is not pinned at w = can be attributed to the strong ex- 
citation of the system with respect to the ground state. A 
quantitative analysis of the spectrum requires the knowl- 
edge of the equilibrium spectrum at the effective temper- 
ature Teff [cf. Eq. dni), Teff = 0.84 for U = 3.3]. Equihb- 
rium spectra are usually computed from imaginary-time 
correlation functions using (maximum entropy) analyti- 
cal continuation, which is not accurate enough at high 
frequencies to allow for a comparison of two rather sim- 
ilar spectra. Using nonequilibrium DMFT, however, we 
can avoid this complication and compute real-time equi- 
librium Green function G^g{t,t') = g^{t — t') without 
analytical continuation (cf. Sec. I VI A|) . Within numer- 
ical accuracy, the resulting equilibrium function indeed 
agrees with the retarded Green function G^{t,t') after 
relaxation (Fig. [7^), which proves that the rapid ther- 
malization at [/ w 3.3 can also be seen in the spectral 
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FIG. 8: Same as Fig. [T] but for the interaction quench to 
U = 5. Spectra (b) are obtained from Fourier transformation 
of real-time quantities, and the Fourier integral (|78|l is cut off 
at Smax = 2.5 with an additional Gaussian factor (see text). 
The corresponding kernel [Eq. (|79|) . k = 0.4] is shown as thin 
solid line. 



function. 

The analysis of the spectrum can now be repeated 
for quenches to the weak- and strong-coupling regime. 
For U <C X^, however, the spectrum remains close to 
the semielliptic density of states for all times, such that 
rather high numerical accuracy would be needed for a 
systematic investigation of the small differences. In the 
strong-coupling regime, on the other hand, the restriction 
to small times t < tmax turns out to be more limiting for 
an investigation of the retarded Green function than for 
static quantities, simply because G^{t + s,t) is known 
only for t < tmax ~ s and not for t < t,nax- Nevertheless, 
one can see that the relaxation of the Green function after 
a quench to U = 5 (Fig. [5^) roughly follows the oscilla- 
tory behavior of the momentum distribution (Fig. \^): 
Close coincidence with the thermal function is reached 
around the time when the jump of the momentum occu- 
pation has its first minimum {t — 0.6), after which the 
deviations to the thermal Green function slightly increase 



again. Similar behavior was found for the double occupa- 
tion, which comes closest to the thermal value at its first 
minimum around t = 0.6;^ In spite of the large effective 
temperature (Tcff = '2V), the spectral function has a clear 
minimum at a; = 0, and well-pronounced Hubbard bands 
at a; w ±11/2 (Fig. [SJd). However, the absolute changes 
with time are small in the strong-coupling regime. This 
behavior is expected because it can be shown that the 
spectrum is independent of time t after a quench to the 
atomic limit. 



D. Optical conductivity 

The two-time optical conductivity a{t, t') describes the 
linear response of the electrical current in a nonequi- 
librium state to a time-dependent electrical field SE{t) 
(which we call the probe field), 

t 

suit)) ^ Jdta{t,t)6E{t). (80) 

— oo 

(Tensor notation of a{t, t') is suppressed.) In solids, opti- 
cal spectroscopy on nonequilibrium states is usually per- 
formed within the pump-probe setup, where the system 
is driven out of equilibrium by a strong laser pulse (the 
pump). In the following we calculate <j{t,t') after the 
interaction quench to sec how the electrical response be- 
comes stationary while the system relaxes towards its 
thermal equilibrium state. 

Microscopically, the optical conductivity is related to 
the current-current correlation function, which can be 
computed from two diagrammatic contributions: (i) The 
bubble diagram of two Green functions Gk and the cur- 
rent vertex Vk = dek/dk, and (ii) diagrams containing 
the vertex corrections of the current vertex.— Within 
equilibrium DMFT, vertex corrections are local and thus 
do not contribute to the conductivity because is anti- 
symmetric under inversion of fc, and Gk is symmetric 
In a nonequilibrium situation these conditions can be vi- 
olated, e.g., due to an electrical pump field, in which 
case the conductivity depends on the relative polariza- 
tion of pump and probe, so that vertex corrections do 
contributCi ^^i^^ However, for the interaction quench the 
inversion symmetry of the state is preserved, and a{t, t') 
can be calculated from the bubble diagram alonci^ 

The microscopic derivation of a(t, t') within nonequi- 
librium DMFT was discussed in detail in Ref. [H. In the 
following we thus only state the results for (7{t, t') after an 
interaction quench in the Hubbard model on the hypercu- 
bic lattice in = oo, with hopping amplitudes that yield 
a semielliptic density of states^ (Eq. ([M]) with V = 1, as 
above). The band dispersion enters the expression via 
the current vertex Vk = de{k)/dk; this is where the hop- 
ping amplitudes enter in addition to the density of states. 
Conductivity is measured in units of (Tq = 2pa^e^V jY? , 
where a is the lattice constant, and p is the number of 
lattice sites per volume. 
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FIG. 9: Optical conductivity a{t + s, t) [Eq. (|80}] after 
quenches to ?/ = 2 (a), !7 = 3.3 (b), and !7 = 5 (c). The 
inset shows a{t,t), and black solid lines correspond to the op- 
tical conductivity in thermal equilibrium at T^a = 0.37 (a), 
Teff = 0.84 (b), and T^s = 2 (c). 



In Fig. [9l cr{t + s,t) is plotted as a function of time- 
difference s. This parametrization is most convenient 
for analyzing how the electrical response of the system 
becomes stationary (i.e., independent of t) during the 
relaxation. The results are compared to the optical con- 
ductivity (Toq(s) in thermal equilibrium, which is ob- 
tained directly from noncquilibrium DMFT without an- 
alytical continuation (cf. Sec. IVI The more familiar 
frequency-dependent optical conductivity 

oo 

cT,^{uj) ^ Rc jdse'^'cT,^{s) (81) 



is plotted in Fig. [TOl 

After quenches to weak-coupling {U = 2, Fig. [5^), 
(j{t,t') undergoes a rapid initial relaxation, but it does 
not approach the thermal value within the accessible 
times. This behavior reflects the prethermalization that 
is observed in the momentum occupation. The conduc- 
tivity at the corresponding effective temperature (Toff = 
0.37) consists of a Drude peak at w = (Fig. fTO]) . which is 
only slightly broadened due to temperature and interac- 
tion. Because a narrow Drude peak implies a slow decay 
of crcq(s) with time difference, we cannot resolve the true 
width of the peak from data which are restricted to small 
times. 

For a quench to U = 3.3 (Fig. Wp)i we observe a rapid 
relaxation of the optical response. The optical conductiv- 
ity depends only on time difference for t^^l/V and coin- 
cides with croq(s) for the effective temperature Toff — 0.67. 
The latter falls off rather quickly with time differences s, 
indicating that the Drude peak is strongly broadened be- 
cause of the large temperature and the relatively strong 
interaction (Fig. fTU]). 

Finally, for the quench to C/ = 5 (Fig. [5J:) relaxation to 
the thermal state becomes again slower than ai U = 3.3. 
We observe the characteristic collapse and revival oscilla- 
tions when a{t + s,t) is plotted at fixed time difference s 
(inset in Fig. [5};). Due to the large effective temperature 
(Tcff = 2) the conductivity of the corresponding equilib- 
rium state is rather a bad metal than an insulator, but 
nevertheless the Hubbard band aX oj = U is clearly sepa- 
rated from the broad feature at = (Fig. [TU|) . 

VII. CONCLUSION 

In this paper we described in detail how weak-coupling 
continuous-time quantum Monte Carlo (QMC) can be 
used as an impurity solver within noncquilibrium DMFT. 
The formalism, which was used in Ref. [l^ to investigate 
the interaction quench in the Hubbard model, was ex- 
tended to the case when the initial state is a finite tem- 
perature equilibrium state at nonzero interaction U . Be- 
cause noncquilibrium experiments in interacting systems 
often start from correlated initial states rather than the 
noninteracting ground state, this extension is a prerequi- 
site to apply DMFT within a variety of experimental sit- 
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FIG. 10: Frequency-dependent optical conductivity in ther- 
mal equilibrium at the temperature T^s = 0.37 {U — 2), 
Toff = 0.84 {U = 3.3), and Tefi = 2 ([/ = 5). The latter 
two curves are scaled by a factor 4 and 8, respectively. All 
curves were obtained by Fourier transformation (|81|l of real- 
time data. In the case oi U = 2, the Fourier integral is cut off 
at Smax = 5 with an additional Gaussian factor, as explained 
for the spectral function. The corresponding kernel [Eq. (|79|l . 
K — 0.2] is shown as dotted line. 

uations in the field of cold atomic gases and timc-rcsolved 
spectroscopy on correlated solids. 

We used the numerically exact QMC solution of the 
DMFT equations to benchmark the generalization of the 
iterated perturbation theory (IPT) to the Keldysh con- 
tour. We find that IPT is remarkably good at weak in- 
teractions. However, in contrast to the equilibrium case 
it yields unphysical results in the intermediate-coupling 
regime and thus cannot provide a reasonable interpola- 
tion between the weak- and strong-coupling regime. The 
reason is that IPT is not a conserving approximation, 
which can lead to an explicit violation of the energy con- 
servation as a function of time in some parameter regime. 

Furthermore, wc used the nonequilibrium formalism 
to solve a system in thermal equilibrium. In this way 
one can avoid analytical continuation and obtain dynam- 
ical quantities in real time instead of imaginary time. 
We used this approach to compute the spectral func- 



tion and the optical conductivity of the single-band Hub- 
bard model. Due to the dynamical sign problem of 
QMC one is restricted to relatively short times, such that 
frequency-dependent quantities, which are obtained from 
real-time functions by Fourier transformation, are consid- 
erably broadened. The real-time formalism can thus not 
directly replace the conventional analytical continuation 
from Matsubara to real frequencies. However, since the 
kernel which mediates the broadening of the spectra is 
explicitly known, it may be useful either to judge the ac- 
curacy of analytically continued spectra, or improve the 
analytical continuation in some frequency range. 

In the last part of this paper we presented further re- 
sults for the interaction quench in the Hubbard model. 
In particular, we investigated the time evolution of the 
real-time Green functions. It was shown that the differ- 
ent relaxation behavior at weak, strong and intermediate 
coupling, which was characterized by the time evolution 
of the double occupation and the momentum distribution 
in Ref. [13, is also reflected in the nonequilibrium spec- 
tral function: In the weak- and strong-coupling regime 
a thermal state cannot be reached within the accessible 
times, whereas the spectrum (as well all quantities that 
can be obtained from it) rapidly relaxes to the thermal 
equilibrium at intermediate coupling (U = Udyn)- 

The fact that the very sensitive f7-dependence of the 
relaxation behavior is manifest also in the spectral func- 
tion suggests that the phenomenon of fast electronic ther- 
malization near Udyn may also be observed with pump- 
probe spectroscopy on correlated systems. Further de- 
tails of this transition-like phenomenon will hopefully 
soon be clarified by means of the DMFT-fQMC formal- 
ism presented in this work. 
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